En el flujo con una cavidad rectangular, el fluido se encuentra confinado en una región rectangular de dimensiones $L\times W$ donde el movimiento se debe al desplazamiento transversal de una de las caras. En este sentido, la velocidad de la cara superior es de $U_0$ y el de las demás interfaces es de $0$. De acuerdo con esto, se tienen las siguientes consideraciones:
De aquí, las ecuaciones de continuidad y de Navier-Stokes se reducen a:
\begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho \vec{u})=0 \quad &\to\quad \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 \\ \rho\left(\frac{\partial \vec{u}}{\partial t}+\vec{u}\cdot\nabla\vec{u}\right)=\rho\vec{f}-\nabla p+\mu\nabla^2\vec{u}\quad &\to\quad \begin{cases} \rho \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) \\ \rho \left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial v}+\mu\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial^2v}{\partial y^2}\right) \end{cases} \end{align}Se puede ver entonces que la ecuación resultante es un caso particular del problema de convección-difusión estable:
$$\frac{\partial f}{\partial t} + \vec{u}\cdot\nabla f=k\nabla^2f+\frac{s}{\rho c_p}$$Donde $\frac{\partial f}{\partial t}=0$, $k = \frac{\mu}{\rho}=\nu$, y $\frac{s}{c_p}=-\frac{\partial p}{\partial *}$ para $f=u, *=x$ y $f=v, *=y$, y las condiciones de frontera corresponden a valores conocidos (tipo Dirichlet).
Ahora, para aplicar el método de elementos espectrales, es necesario realizar los siguientes pasos:
Cada uno de los pasos se tratará en detalle a continuación
In [1]:
# Prepara el ambiente de trabajo con las librerías a utilizar
import numpy as np # Ofrece funciones para operaciones básicas con arreglos
import scipy.linalg as lin # Ofrece funciones para operaciones de álgebra lineal
import matplotlib.pyplot as plt # Permite incluir un ambiente de visualización
from matplotlib.patches import Patch
from matplotlib.collections import PatchCollection
from mpl_toolkits.mplot3d.axes3d import Axes3D
import timeit # Permite la medición de tiempos para los distintos algoritmos
from sem2D import * # Agrupa las funciones externas en un archivo
# from __future__ import division # Corrige posibles problemas que surgen a la hora de dividir enteros
#%pylab inline
In [2]:
## Valores de prueba
L = 2.4
W = 1.6
mu = 1.
U0 = 1.
In [3]:
Ne = np.array([12,8]) # Indica la cantidad de elementos a lo largo de los ejes x y y
m = 3 # Indica la cantidad de divisiones sobre cada elemento en el plano xi-eta (admite valores 1, 2 y 3)
# Genera los puntos correspondientes a los vertices de cada rectángulo
(px, hx) = np.linspace(-L/2,L/2,Ne[0]+1,retstep=True)
(py, hy) = np.linspace(-W/2,W/2,Ne[1]+1,retstep=True)
Conociendo el número de elementos y su distribución, es posible generar los nodos correspondientes a los vértices de cada elemento. Sin embargo, como no se trabaja con elementos lineales, es necesario generar los nodos de interpolación para cada elemento. Para ello, es necesario definir la familia de polinomios sobre la cual van a estar basados estos nodos. En este caso se utilizan los nodos de interpolación de Lobatto.
In [4]:
(x, y, ne, npe) = genNodesTri(Ne,m,px,py)
# Grafica la malla de puntos
for l in range(ne):
plt.plot([x[l,0],x[l,m],x[l,2*m],x[l,0]],
[y[l,0],y[l,m],y[l,2*m],y[l,0]],'r--')
plt.hold(True)
plt.plot(x,y,'bo',ms=4)
plt.xlabel('coordenada x')
plt.ylabel('coordenada y')
plt.show()
In [5]:
h = min(hx,hy)
(coords, C, gfl, ng) = genConn(Ne,m,ne,npe,x,y,L,W,U0,h,tipo='tri')
Para la generación de la matriz del problema es necesario definir las matrices de difusión, masa y advección por cada elemento:
$$A_{ij}^{(\ell)}=\iint_{E_\ell} \nabla\psi_i^{(\ell)}\cdot\nabla\psi_j^{(\ell)}dxdy \quad B_{ij}^{(\ell)}=\iint_{E_\ell} \psi_i^{(\ell)}\psi_j^{(\ell)}dxdy \quad C_{ij}^{(\ell)}=\iint_{E_\ell} \psi_i^{(\ell)}\vec{u}\cdot\nabla\psi_j^{(\ell)}dxdy$$Estas integrales se obtienen utilizando la cuadratura de Gauss triangular aprovechando la geometría sencilla del dominio. Dado que el orden máximo de las funciones de interpolación es 3, bastaría con usar la cuadratura de orden 4. Puntualmente, usando las funciones ortogonales de Proriol las funciones de interpolación se pueden describir como:
$$ \psi_i(\xi,\eta) = \frac{\text{Det}[V_\phi(\xi_1,\eta_1,\xi_2,\eta_2,\dots,\xi,\eta,\dots,\xi_N,\eta_N)]}{\text{Det}[V_\phi(\xi_1,\eta_1,\xi_2,\eta_2,\dots,\xi_i,\eta_i,\dots,\xi_N,\eta_N)]} $$Comenzamos con la matriz que define el sistema de ecuaciones, o la matriz de difusión. Dado que esta depende de los gradientes de las funciones de interpolación, es necesario realizar todos los pasos de construcción de los términos como sigue.
De acuerdo con esta descripción, se tiene entonces que las derivadas con respecto a $\xi$ y $\eta$ se obtienen derivando los polinomios de proriol. Esto radica en que las únicas entradas no constantes son las que corresponden a las entradas $\xi$ y $\eta$ del numerador.
De este modo, la función de transformación de coordenadas corresponde a:
$$ \vec{x} = \sum_{i=0}^{N_p-1} \vec{x}_i^E \psi_i(\xi,\eta) = \sum_{i=0}^{N_p-1} \vec{x}_i^E \psi_{(j,k)}(\xi,\eta) = \sum_{i=0}^{N_p-1} \vec{x}_i^E \mathcal{L}_j(\xi)\mathcal{W}_k(\eta) $$y por lo tanto sus derivadas son:
\begin{align} \frac{\partial x}{\partial \xi} &= \sum_{i=0}^{N_p-1}x_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &\quad \frac{\partial x}{\partial \eta} &= \sum_{i=0}^{N_p-1}x_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} \\ \frac{\partial y}{\partial \xi} &= \sum_{i=0}^{N_p-1}y_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &\quad \frac{\partial y}{\partial \eta} &= \sum_{i=0}^{N_p-1}y_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} \end{align}Ya con estos valores es posible soluciones los gradientes de las funciones solucionando los sistemas:
$$ \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{pmatrix} \nabla\psi_i = {\frac{\partial\psi_i}{\partial\xi} \choose \frac{\partial\psi_i}{\partial\xi}} $$De manera que el coeficiente de cambio de coordenadas corresponde a:
$$ Det[J] = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{vmatrix} = \frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta}-\frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi} $$Ahora bien, como tampoco se conocen los valores de presión, es necesario incorporarlas también como incógnitas. Para ello se calculan los vectores: $$ D_{i}^{(\ell)}=\iint_{E_\ell} \nabla\psi_i^{(\ell)}dxdy$$
Para el otro lado de la ecuación se implementa la función de entrada. Para esto se recorre sobre todos los elementos y se van calculando los valores para la función de manera correspondiente. Estos valores se acomodan al sistema global aprovechando la matriz de conectividad.
Así mismo, es necesario modificar la matriz global de manera que se incorporen las condiciones de frontera de Dirichlet. Si $l$ es un nodo de frontera, esto se hace a través de los pasos:
Así, se utilizan las funciones auxiliares (precalculando los valores para los polinomios requeridos) siguientes:
In [6]:
(Mat_dif, vec_b) = matDiffTri(Ne,m,ne,npe,ng,x,y,C,gfl,mu)
plt.spy(Mat_dif)
plt.show()
Teniendo ya construidos los componentes del sistema a resolver, basta aplicar un esquema de solución razonable para el tipo de matriz obtenida. En este caso se utiliza el solucionador por defecto del módulo linalg. Dado que los valores para el primer y último nodo son conocidos, se omiten los valores correspondientes.
In [7]:
start = timeit.default_timer()
sols = lin.solve(Mat_dif[1:,1:],vec_b[1:])
print timeit.default_timer()-start
# Recupera las velocidades
solx = np.concatenate([[gfl[0,1]], sols[0:ng-1]])
soly = sols[ng-1:2*ng-1]
# Calcula la magnitud de la velocidad en cada nodo
solm = np.abs(solx**2+soly**2)
# Recupera las presiones
solp = sols[2*ng-1:]
In [8]:
# Grafica las fronteras de los elementos
for l in range(ne):
plt.plot([x[l,0],x[l,m],x[l,2*m],x[l,0]],
[y[l,0],y[l,m],y[l,2*m],y[l,0]],
'r--')
plt.hold(True)
# Grafica el campo de velocidades
plt.quiver(coords[:,0],coords[:,1],solx,soly,color='k')
plt.title('Campo de Velocidades para el Flujo en la Cavidad')
plt.xlabel('Posicion (coord x)')
plt.ylabel('Posicion (coord y)')
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])
plt.show()
In [9]:
# Grafica el campo de presiones
xpres = np.zeros(ne)
ypres = np.zeros(ne)
for l in range(ne):
xpres[l] = (x[l,0]+x[l,m]+x[l,2*m])/3.0
ypres[l] = (y[l,0]+y[l,m]+y[l,2*m])/3.0
zpres = solp-np.mean(solp) # Como se trabaja con el gradiente es importante
# la manera en que cambia, no su valor puntual
fig = plt.figure()
ax = fig.gca(projection='3d')
# Grafica las fronteras de los elementos
for l in range(ne):
plt.plot([x[l,0],x[l,m],x[l,2*m],x[l,0]],
[y[l,0],y[l,m],y[l,2*m],y[l,0]],
'r--')
ax.plot([xpres[l],xpres[l]],[ypres[l],ypres[l]],[0,zpres[l]],'b',linewidth=2)
plt.hold(True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Presion')
plt.show()
In [ ]: